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Stretched polymers with attractive interaction are studied in two and three dimensions. They 
are described by biased self-avoiding random walks with nearest neighbour attraction. The bias 
corresponds to opposite forces applied to the first and last monomers. We show that both in d = 2 
and d = 3 a phase transition occurs as this force is increased beyond a critical value, where the 
polymer changes from a collapsed globule to a stretched configuration. This transition is second 
order in d = 2 and first order in d = 3. For d = 2 we predict the transition point quantitatively 
from properties of the unstretched polymer. This is not possible in d = 3, but even there we can 
estimate the transition point precisely, and we can study the scaling at temperatures slightly below 
the collapse temperature of the unstretched polymer. We find very large finite size corrections which 
would make very difficult the estimate of the transition point from straightforward simulations. 



I. INTRODUCTION 

Deformations of polymers have been studied for many 
years, such as stretching a single DNA or titin macro- 
molecule • Such experiments will become even more 
important with the rapid improvement in single molecule 
experiments |0| . Their theoretical understanding has at- 
tracted much attention. Most such experiments have 
been performed in good solvents, where also the theory is 
best understood [|| . But of particular interest in biology 
is the stretching of a collapsed polymer, i.e. of a polymer 
in a poor solvent below the 0-temperature. Up to now 
there seem to exist only few experiments in this regime 
1^. Unfolding proteins in this way could e.g. give impor- 
tant information on their spontaneous folding pathways. 

Although there exist several theoretical papers on 
stretching of collapsed polymers 0-^, we beheve that 
further work is needed for a full understanding. It is 
generally believed that there is a first order transition 
between the globule and the stretched phases in 3 dimen- 
sions. But it seems that precise estimates of the critical 
force do not exist. Also, there are no predictions for the 
scaling laws expected when the temperature approaches 
the collapse temperature Tg from below. Finally, it seems 
that there is only one paper which deals with d = 2. 
Moreover, while it is claimed in that the transition 
is also first order in d = 2, we shall find there a rather 
different situation and a second order transition. 

A stretched polymer in a bad solvent is modeled as 
a biased interacting self-avoiding random walk (BISAW) 
on a regular lattice (square in d = 2, simple cubic in 
d = 3) with nearest neighbour attraction. In this model 
no two monomers can visit the same site. The attraction 
is taken into account by a Boltzmann factor q = e~^'^ 
for each pair of non-bonded monomers occupying near- 
est neighbours on the lattice. Here /3 — 1/kT and e < 
is the attractive potential between non-bonded nearest- 
neighbour pairs. The stretching is described by a factor 
b = exp{f3aF) where a is the lattice constant and F is 
the stretching force Q . The partition sum is therefore 



z=Y. rb\ (1) 

walks 

where m is the number of non-bonded occupied nearest- 
neighbour pairs, and x is the distance (in units of lattice 
constants) between the two end points of the chain in the 
direction of F. We are interested in polymers which form 
collapsed globules if they are not stretched, i.e. q > qg 
where qg = e"'^/*^^" . 

In the present paper we employ the Pruned-Enrichcd- 
Rosenbluth Method (PERM) to study the phase 
transition of the BISAW in d = 2 and 3. But we shall 
see that, in order to understand the BISAW, we need 
also some other results. More precisely, we will need 
also some more results on the unbiased interacting self- 
avoiding random walk (ISAW), both on infinite and on 
finite lattices. On finite lattices with periodic boundary 
conditions we can study the bulk behaviour of collapsed 
polymeric matter without being disturbed by surface ef- 
fects |9|,[8). 

As we shall see, a stretched collapsed polymer in d = 2 
forms, in the infinite chain length limit TV — > cxd and for 
forces below the transition to the stretched phase, a com- 
pact object of a shape shown in Fig. 1. There the bound- 
aries are circular arcs, the density inside is independent 
of the stretching force F, and the shape (i.e. the angles 
at the two extreme points) is determined both by F and 
by q. The transition to the stretched phase occurs when 
these angles tend to zero. 

In d = 3 the situation is more complicated because of 
the Rayleigh instability. For F < Fc one has an object 
of roughly elongated ellipsoidal shape, but for F — Fc 
two phases coexist: Part of the chain is stretched, while 
the rest is still collapsed (Fig. 2). In the collapsed part, 
the density is the same as in an unstretched globule (for 
N oo). The stretched part might consist of a sin- 
gle piece (tadpole configuration pO[), or of two pieces as 
in Fig. 2. Within our theory both configurations have 
the same energies. The critical force Fc vanishes linearly 
when T Tg. 
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The paper is organized as follows: In Sec. II we treat 
the 2-dimensional case. We first formulate in more detail 
the model sketched above, and we show that it is in good 
agreement with simulations both of stretched and un- 
stretched chains. The 3-dimensional case are discussed 
in Sec. III. Conclusions are drawn in Sec. IV. Finally, 
details of the Monte Carlo algorithm are given in an ap- 
pendix. 
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FIG. 1. (a) The geometric shape of a stretched polymer in 
a poor solvent in d = 2, (b) Coordinates used in the calcula- 
tion. 




FIG. 2. Schematic drawing of a stretched polymer in a poor 
solvent in d = 3, at an intermediate force F where the col- 
lapsed and extended phases coexist. 



II. D=2 

Let us first consider an unstretched polymer in d di- 
mensions, for any d > 2, below the 0-point. In the follow- 
ing, the lattice constant will always be a 1. For large 
N and for infinite lattices the polymer forms a globule 
with a bulk monomer density p{T) which tends to zero 
for T Te, p{T) - [Te - Tf. For d > 3 we have {3 = 1 
and logarithmic corrections , but for d = 2 the power 
/3 seems to be unknown. For T ^ we have p{T) 1. 



The free energy for such a polymer consists of two 
parts: An extensive bulk contribution cx N, and a surface 
contribution cx iV^''^^^/'^. For d = 2 this reads 

-\0gZN{q,b = l)^ Po.{q)N + a{q)^ (2) 

with Poo being the chemical potential per monomer in 
an infinite chain, and a is related to the surface tension 
(=free energy per unit of perimeter length) a by 

ay/N = 2ttR(j (3) 

with ttB? p = N . Combining the last two equations we 
obtain 

Now assume we start stretching the polymer by pulling 
at the two ends. The first effect will be that the ends are 
drawn toward the surface and the whole globule is ro- 
tated such as to minimize the potential energy. We will 
neglect the change in entropy associated to that. When 
the stretching force increases, we develop a shape as in- 
dicated in Fig. la. We shall assume that the bulk is 
incompressible, so that we still have the same bulk free 
energy, and we obtain our final ansatz 

~\ogZN{qM^ l^ooN + (jP-2LF (5) 

where P is the perimeter and 2L is the end-to-end dis- 
tance. Parameterizing the shape by a function y = y{x) 
with ?/(0) = yo, dy/dx\x=o = 0, and y{L) = (see 
Fig. lb), we have 

P = 4: dxy^l + {dy/dxY, 4p / dxy{x) = N. 
Jo Jo 

(6) 

The optimal shape is obtained by maximizing Zn at 
constant N . The Euler equations for this maximization 
problem consist of two parts. The first, resulting from 
dZN{q,b)/dL — 0, gives the equilibrium condition be- 
tween stretching and surface tension forces, 

F = 2(jcos(j), (7) 

where the angle (f) is indicated in Fig. la. Thus the criti- 
cal point, corresponding to = 0, is at — 2cr or 

be = e^- . (8) 

The second part, maximization with respect to y{x), 
gives that y{x) is a circular arc. Using the angle indi- 
cated in Fig. la, we have then 

P^4L^, N/p^2L^^^^^4^. (9) 
sin sin 

For b<bcwe have 
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log ZN{q,b) = /XooiV + (J^^/^[8(0- sin0cos(/))/p]i/2 

(10) 



with 



arccos ■ 



log 6 



(11) 



Eqs. form our final solution. They involve the 

three temperature dependent material constants /Xqo, ct, 
and p. The former two can be estimated from Eq. (||) 
if we measure in addition the gyration radius, but more 
precise estimates of /Zoo and p result from simulations on 
finite lattices in the dense limit ||2l|,|l^. Thus our strat- 
egy in verifying the above theory numerically consists of 
the following steps: 

1. We simulate chains on finite lattices of size L x L 
(typically with i = 8 to 64) by means of the PERM 
algorithm which gives us directly estimates of the 
partition sum. For q > qg its logarithm is not 
convex. For suitably adjusted pL, log2'jv(q, 6 = 
1,L) + plN has two peaks of equal heights: One 
at « 0, the other at iV « L^p. Extrapolating to 
L ~ oo we obtain pao- 

2. In the next step we simulate chains on (practically) 
infinite lattices, again by means of PERM. Using 
Eq. (H) and the already obtained values of poo and 
p we then obtain a. 

3. Finally we simulate stretched polymers, to com- 
pare with the prediction of Eq. (p^. In addition to 
Zjf(q,b) we measure in these runs also (x) which 
should be equal to 2L for large N as long as 6 < be- 
We should, however, immediately warn that these 
latter measurements are not very conclusive, since 
finite size corrections are large in this regime. We 
might point out that for 6=1 there is the exact 
result 



(12) 



Since (x^) ^ N for collapsed polymer in d = 2, we 
thus have (x) ~ N for very small b. This shows that 
the above model (which would give (x) ~ 
cannot be correct for small N, as we had indeed 
pointed out already before. 

In Fig. 3 we show log + plN for q = 2.4, b = 1, and 
for finite lattices with L"^ = 2"^, 2*, . . . , 2^^ sites, plot- 
ted against N/L"^. The curves are in the same order as 
shown in the legend. Since qg « 1.95 for this model [ p2[ , 
this is deep inside the collapsed regime. We used helical 
boundary conditions. When was not an integer, the 
effective lattice shape was not a perfect square, but this 
led to negligible corrections. The values oi plN were 



fixed by demanding the right hand peak to have height 
exactly zero. The values obtained this way are shown in 
Fig. 4. They are plotted there against i^^-^-^ because 
this gave the best straight extrapolation to L = cxd. The 
extrapolated value is poo — —1.3213(1). The positions of 
the right hand peaks give finite lattice approximations to 
the density p. Extrapolating to L = oo gave p — 0.84(1). 
Analogous simulations were also done at different values 
of q. 
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FIG. 3. For d = 2, g = 2.4 and fo = 1, In Zjv -f plN versus 
N/L'^ for finite lattices with L'^ = 2^ , 2*, . . . , 2^^. The values 



ol PL were fixed by demanding that the peaks at N/L 
have zero height. 
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FIG. 4. PL as obtained from Fig. 3, plotted against L~^ ''^ 
which gave the most straight extrapolation. Error bars are ca. 
0.00005, smaller than the symbol sizes. The extrapolation to 
i = oo is ^oo = -1.3213 ± 0.0001. 



Results from unbiased ISAW on an 'infinite' lattice 
(i.e. on a lattice that was so large that no walk reached 
the boundary) are shown in Fig. 5. There we plot 
log^AT -I- pooN against N^^"^. From Eq. (^) we expect 
this to give a straight line with slope —a. Actually the 
line is slightly curved, indicating that there are further fi- 
nite size corrections for small N and systematic sampling 
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corrections at large N. The latter are indeed to be ex- 
pected. When we constructed histograms of tour weights 
(see the appendix), we found that the simulations are un- 
reliable for > 1000. We should have shown only the 
data for N < 1000, but we showed all data, for the fol- 
lowing reasons: 

• we want to compare these simulations with biased 
and finite volume simulations of the same length 
which are easier for PERM and which are thus still 
reliable for N = 3500, the longest chains used in 
Fig. 5. 

• Even if we cannot be sure that the data for N > 

1000 arc correct, wc cannot argue cither that they 
must be wrong. They are most likely too low, 
since PERM has difficulties to sample configura- 
tions whic;li start out (as the chain grows) to have 
high energy (i.e. few contacts), but which become 
'good' during later growth stages. But we do not 
know how important this is. 

The curvature observed in Fig. 5 dominates the error 

in our estimate a = 0.46(2). Combining this with the 
previous estimate for p we obtain a = 0.119(5), and from 
this we predict be = 1.269(13). Remember that all this 
is for q = 2.4. 

Finally, results from simulations with b > 1 are shown 
in Figs. 6 to 9. In Fig. 6 we show log ZN{q, b) + Hoo{<l)N 
versus N. We expect these curves to become horizontal 
for ^ oo as long as 6 < 5c, while they should increase 
linearly for b > be- This is indeed seen, although we now 
find a slightly larger value be = 1.285(6). In spite of the 
small discrepancy with the predicted value we consider 
this as remarkable agreement. Our data are not precise 
enough to allow a comparison with the detailed predic- 
tions for 1 < 6 < 6c. 

For b > be the chemical potentials n'^^{q,b) in the 
stretched phase can be estimated easily by demanding 
that logZjv(Q', 6) -|- iJ,'^{q,b)N becomes Af- independent 
for large N. The results are shown in Fig. 7. Although 
this figure clearly shows a continuous transition, the data 
are not precise enough to quote a meaningful critical ex- 
ponent. It is compatible with 3/2, but the uncertainty is 
large. 

Among the other measurements during these runs, the 
most interesting are those of (x). Again they are not 
precise enough for a detailed comparison with the pre- 
dictions in the regime b < be- But they also show clearly 
the phase transition at 6 = be, since; (.;•) rises linearly 
with N only in the stretched phase (sec Fig. 8). Values 
oi V = d{x)/dN, obtained by extrapolating the observed 
slopes in plots of (x) against N toward N = oo, are shown 
in Fig. 9. Again we see a continuous transition with an 
exponent which is roughly equal to 0.7, but uncertainties 
are too large to make a more definite statement. 




FIG. 5. Values of \nZN(q,b = 1) + Moo^V versus N^^^ for 
I = 2, q = 2.4, and /Xtxj = —1.3213. The dashed line has slope 
-a (ct = 0.46 ±0.02). 
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FIG. 6. Values of In ZN{q, b) +ij,ooN versus A'' for the same 
q and //oo as in Fig. 5, and for various values of 6. 
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FIG. 7. Values of \imN^^[N-'^ In ZN(q,b)] versus &, for 



d = 2 
b < be 



and q = 2.4. For b > be this is —n'£{q,b). For 
(~ 1.3), wo replaced the actual numerical estimates 



by the theoretical value — Atoo(g). 
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FIG. 8. Average displacement (x) in the bias direction for 
= 2, q = 2.4 and for various values of b. 
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FIG. 9. Asymptotic "velocity" d{x)/dN in the bias direc- 
tion for d = 2, q = 2.4 and for various values of b. 



Before leaving this section we should remark that we 
checked carefully, in view of Ref. jl^, that the unfold- 
ing transition is indeed second order in our model. We 
looked e.g. at histograms analogous to those in Figs. 12 
and 13. We definitely saw no hint of bistability or hys- 
teresis. Together with the good agreement with the above 
model this seems definitely to rule out a first order tran- 
sition. We cannot rule out, however, that the transition 
is first order in another microscopic model such as the 
bond fluctuation model studied in 1161 . 



III. D = 3 

Some of the treatment for d = 2 can be carried over 
to 3 dimensions with minor modifications. Instead of 
Eq. (|) we now have ||2|,|l|| 

- logZNiq, - 1) « f^oo{q)N + aiq)N^^^ . 



(13) 

Again we estimate fioc most reliably from finite lattice 
simulations in the dense limit [n8[ . Again a can be related 



to a surface tension, and can be estimated numerically 
by plotting logZjv(g, 6 = 1) -|- ^aoiq)N against N'^^^. 

In the weak stretching regime we could still try to solve 
the optimal shape, but this time the minimization prob- 
lem is more complicated and we were not able to find 
an explicit solution. But more important, we do not ex- 
pect a continuous transition as in d = 2. The reason is 
that we expect a Rayleigh instability when the globule is 
stretched too much, which is expected to occur, before 
the limit (/< = is reached. Thus we cannot expect to be 
able to predict the transition point as in d = 2. 

On the other hand, things also simplify because the 
critical exponents at the Q point now should be mean 
field like (the upper critical dimension for the Q point is 
3 [^). This means in particular that (T,p and ^oo{q) — 
Mcio(<Zc) all should vanish as ^ {Tg — T) when T ^ Tg 
from below. 

For a stretched globule with & < 6c we expect that 
Zn differs from the value for b — 1 only by its changed 
surface, i.e. by terms ^ N^^^. For b > be we have a 
string-like phase whose free energy should increase lin- 
early with N with very small corrections, 



logZNiq,b>be)^ ^i<^'Hq,b)N , 



(14) 



where the chemical potential /i*^*-* {q, b) in the stretched 
phase is a function independent of the chemical potential 
/^oo('Z, ^) in the globular phase. The end-to-end distance 
in the stretched phase is given by {x) = Nv{q, b) with 



v{q,b) = -b 



9/i(")((7,6) 



db 



(15) 



Stability of the stretched phase requires that 
dv{q,b)/db > 0, i.e. 



dfi^'''>{q,b) .,9 



db 



(6— )V*)(q,6) <0. 



db 



(16) 



In between these two regimes we expect a coexistence 
region where part of the chain forms a (single) globule, 
while the rest forms one or two stretched pieces. A some- 
what crude model of a polymer in the coexistence region 
which, however, catches all essential features including 
finite size effects, is the following: The stretched part 
has Ns monomers, and its partition sum is described by 
Eq. (|lj) with N replaced simply by N^. The globular 
part has Ng = N — monomers, and its partition sum 
is described by Eq. ( p^ with N replaced by Ng . The total 
partition sum is just the product of the two, i.e. the free 
energy is just the sum of free energies of the two parts. 
Notice that this involves a number of approximations: 

• There is no penalty for the area where the glob- 
ule and the stretched part(s) are attached to each 
other. Such a free energy contribution should be 
independent of N and can be safely neglected. 
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• The globule is approximated by a sphere. This 
is a more serious approximation. It systemati- 
cally overestimates the free energy, by an amount 



c>c Ng^^. This is of the same order of magnitude as 
if a wrong surface tension were used. It should 
therefore lead to quantitative errors, but not to 
qualitative ones. Moreover, the errors should be 
small because too elongated globules are prevented 
by the Rayleigh instability. 

• We neglect all fluctuations. 

The total free energy in the coexistence region is thus, 
for a fixed total end-to-end distance x (and not denoting 
explicitly the dependence on q), 

- log Zn = Ng^ioo + N^'^d + (iV - 7Vg)/i(^) (6) - xF . 

(17) 

Minimizing this with respect to Ng we obtain 

^ioo-^i^'\b)+'^N-^/^ ^Q. (18) 



This can be read as an implicit equation for h. If solved, 
it gives 6 as a function of Ng. Notice that h is indepen- 
dent of Ns : According to our model the force needed to 
pull out the chain from the globule depends on its size 
through the surface tension (last term in Eq. (|l8|)), but 
not on the length Ns of the stretched part. 

In the limit Ng oo this term does not contribute, 
and we obtain the condition for the true transition point 

'^'Hbc) ^ f^oo ■ (19) 



For finite Ng we see that at the coexistence point /z'*-'(6) 
decreases with Ng. Since d fi'^^^ {b) / db < 0, this implies 
that the effective critical b increases with Ng. In other 
words, as the chain is pulled out from the globule, Ng 
shrinks, and thus the force needed to pull out more of 
the chain decreases. This is the basic instability which 
makes the transition first order and shows that our ansatz 
is consistent. It implies that a finite globule will be en- 
tirely pulled open as soon as 6 > b^^ — const/ N^^^ . Notice 
that a more realistic model (where the spherical glob- 
ule is replaced by some other shape with limited aspect 
ratio) would still give the same qualitative results: For in- 
finitely large globules, the critical point is given by equat- 
ing fJ,^^^ {be) with the chemical potential for unstretched 
globules; and finite size corrections to this are negative 
and decrease as N~^^^. 

In a first set of simulations we determined numeri- 
cally ^00(9) and p{q) for a wide range of q. We did 
this again, as in d = 2, by performing the simulations 
on finite lattices in the dense limit. In a second step, 
we obtained a from simulations of unbiased ISAWs in 
infinite lattices. In contrast to d = 2, here we encoun- 
tered the problem that there are significant further cor- 
rections to the asymptotic ansatz Eq. (13). We illustrate 



this in Fig. 10 where we plot log[.ZAr_„((7)/Zjv+„(g)]/(2n) 
against N~^^^ for b = 1 and several values of q. Here, 
n — 1 + [A^/20J. The points on the y-axis are obtained 
from dense limit simulations. The straight lines are ex- 
trapolations constrained to pass through these points. 
Their slopes are 2a /3. While the curves are compati- 
ble with these extrapolations, close inspection shows that 
they are not convex as one might have guessed naively. 
Thus, extrapolations not aided by the dense limit simu- 
lations would be prone to large errors. 
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10. ln[ZN-n{q)/ZN+niq)]/i2n) versus N'^^^ for 
6=1 and several values of q. The points on the 
y-axis are obtained from dense limit simulations on finite lat- 
tices, and the straight lines are extrapolations to ^ 00, 
constrained to pass through these points. 



Finally, we performed simulations of BISAW (on 'infi- 
nite' lattices). We shall discuss in detail only the most 
extensive simulations, done at q = 1.5 which is deep in 
the collapsed region (ge — 1-3087(3) for this model [psf). 
But similar simulations were also done at different values 
of q. 
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FIG. 11. Average displacement in the bias direction (x) for 
d = 3, q = 1.5 and for various values of b. 
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Average displacements (x) against N for various bi- 
ases b are shown in Fig. 11. For b = 1, we have of course 
{x) = 0. With increasing b, (x) still remains close to zero 
for large N and increases very slowly as TV ^ oo as long 
as b is close to 1. But for small N it increases roughly 
oc N , with an abrupt decrease in a narrow range which 
shifts toward larger N when b is increased. As b increases 
above a certain value be, the linear range of (x) versus 
N extends to = cx), indicating the phase transition 
from the collapsed to the stretched polymer phase. From 
Fig. 11 the critical point of the finite system might be 
guessed to be between 1.60 and 1.65, but we will see that 
finite-size corrections are very large and be is actually 
larger than 1.85. 

Very similar behaviour is seen for i?^ and for the aver- 
age number of non-bonded nearest-neighbour pairs (m) 
(not shown). For b > be they increase linearly with N, 
Rj^f cx (m) cx N. For b < be but close to be they first 
show the same behaviour, indicating that short polymers 
would be stretched at this b, but then cross over to the 
collapsed phase in a very narrow range of N. 

We claim that this is a first (yet inconclusive) indica- 
tion for the first order nature of the transition. More di- 
rect indications are obtained by studying the histograms 
of X and m. The histograms shown in the following fig- 
ures are obtained by combining data from runs with dif- 
ferent q and b and reweighting them. Combining MC 
results from different temperatures is not trivial for con- 
ventional Metropolis- type Monte Carlo algorithms where 
absolute normalization is unknown p^ . In contrast, it 
is straightforward for PERM, since PERM gives directly 
estimates of the partition sum and of the properly nor- 
malized histograms 



'5x 



(20) 



walks 



Reweighting histograms obtained with runs performed 
nominally at q and 6 is trivially done by 



P,:.y{m,x) = Pg,,{m,x){q'/qy^{b'/bY 



(21) 



Combining results from different runs can then be either 
done by selecting for each (to, x) just the run which pro- 
duced the least noisy data (which was done here in most 
cases), or by assuming that the statistical weights of dif- 
ferent runs are proportional to the number of 'tours' p^ ] 
which contributed to Pq.b{m, x). 

In the present work we studied only single variable his- 
tograms Pq(rn) and Pb{x), for which the above holds with 
the appropriate modifications. 

Histograms of m for fixed q and b, and for three differ- 
ent values of A, are shown in Fig. 12. For N — 1500 and 
A = 2000 we see two peaks, corresponding to the col- 
lapsed (right hand, m k, 0.9A) and stretched (left hand, 
TO ~ 0.4A) phases. For A > 2000 all chains would be 
collapsed, while for A < 1000 all chains are stretched, in 
agreement with Fig. 11. 
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0.007 




0.6 0.8 1 1.2 1.4 

m/N 

FIG. 12. Histograms of the number of non-bonded near- 
est-neighbour pairs P{m) versus m/N for d = 3, q = 1.5, and 
b — 1.55 < be- The peaks near m/N « 0.4 correspond to the 
stretched phase, the ones near m/N « 0.9 to the collapsed 
phase. Notice that chains with TV = 1000 are for these values 
of b and q entirely in the stretched phase (there is no peak 
near m/N « 0.9), in agreement with Fig. 11. Normalization 
is arbitrary. 
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FIG. 13. Histograms of the end point distance P{x) versus 
x/N for q = 1.5. Biases were adjusted so that both peaks 
have equal height: b = 1.4040 {N = 500), 1.4925 (N = 1000), 
1.5386 {N = 1500), 1.5658 {N = 2000), 1.5855 {N = 2500). 
Normalization is arbitrary. The peak at a;/A'^ « corresponds 
to the collapsed phase, the other to the stretched phase. 



Analogous histograms of x are shown in Fig. 13 for a 
wider range of A. We now see an even more pronounced 
double peak structure, with the left (right) peaks cor- 
responding to the collapsed (stretched) phase. In this 
figure we kept q = 1.5 fixed but varied 6, so that both 
peaks have the same height for each A. In addition we 
adjusted the normalization arbitrarily such that all peaks 
have similar heights. We see clearly that the height of the 
minimum between the peaks shrinks to zero for A — > oo, 
and that the horizontal distance between the peaks in- 
creases with A. Taken together, they form a clear indica- 
tion for a first order transition. Notice that a double peak 
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structure with decreasing minimum alone would not be a 
conclusive proof, as shown by the 0-point in dimensions 
d> 4 [||. 

The values b = bc{N) for which the two peaks to have 
equal height (indicated in the caption of Fig. 13) are ef- 
fective finite N transition points. According to our phe- 
nomenological model we expect them to scale as 

I 6c(iV) - |oc . (22) 

The values obtained from Fig. 13 are plotted in Fig. 14 
against N~^^^. From Eq. (^2|) we expect them to fall 
onto a straight line. This is indeed the case, and the 
extrapolation to — > oo gives be — 1.84 ± 0.04. 



to the transition point, but also much lower entropy, so 
that it can easily be missed during a run with finite CPU 
time. Whenever this happened, were plotted in Fig. 15 
the value predicted by Eq. (|2^) instead of using the direct 
estimate of ZN{q,b). 

Since the curve for 6=1 becomes horizontal for 
N oo, the true phase transition occurs at that value of 
b for which the straight line in Fig. 15 is also horizontal. 
This can be estimated very easily and with high preci- 
sion, giving for q — 1.5 our final estimate be ~ 1.856(1). 
This is in perfect agreement with the above finite- ex- 
trapolation. Its error is dominated by the uncertainty of 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 

1^-0.33 

FIG. 14. EfFective transition points bc{N) versus 

for 

d = 3, q = 1.5. The true transition point be — 1.84 ± 0.04 is 
determined by extrapolating to the y-axis. 
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FIG. 15. In ^jv(q, b) + fiooN for d = 3,q — 1.5, and for var- 
ious values of b. The value /loo = —1.7530 ± 0.0003 used in 
this plot was obtained from dense limit simulations on finite 
lattices. 



This estimate of be is already more precise than any es- 
timate we could obtain from Fig. 11. But an even more 
precise estimate is obtained by analyzing the partition 
sum itself. In Fig. 15 we show log Z^iq, b) -f fiooN versus 
A. For small b the curves are close to the curve ioi b — 1. 
As b increases, the initial (small- A^) parts of these curves 
are straight lines with less and less negative slopes. In 
this regime the polymer is stretched. As long as these 
slopes are negative, the straight lines will intersect the 
curve for 6 = 1 at some finite value of A, say Nc{b). Ob- 
viously these arc the chain lengths where b = bc{Nc) (see 
Fig. 14). 

For A^ > Nc{b), the values of log Zpf{q, b) + fiooN must 
deviate from the straight lines. Indeed, the curves in 
Fig. 15 cannot cross each other since |26) 

log ZN{q,b) - log ZN{q,b = 1) 

= log[(cosh(a;log6))f,=i] 
(log&)2 



> log[l 



(23) 



In the simulations, such a crossing can of course happen 
due to metastability of the stretched phase. With our al- 
gorithm, the collapsed state has much lower energy close 



The results of log(6c('?)) and o-(g) for a wide range of q 
values are shown in Fig. 16. We see that both curves are 
roughly linear near the O point, with slopes close to one. 
Exact scaling laws log{bc{q)) ^ q — qe and a{q) ^ q — qg 
cannot be expected because of the strong logarithmic cor- 
rections at the O point |ig|,[lq|. 




0.01 

0.01 0.1 

(q - q9)/qe 

FIG. 16. Log-log plots of lnfoc(<?) and of surface tensions 
a(g) versus (g — qe)/qe for d = 3. The continuous lines are 
drawn to guide the eye. 
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IV. CONCLUSION 

By employing the PERM algorithm to the BISAW 
model with attractive interaction on square and simple 
cubic lattices, we have studied the process of stretching 
collapsed polymers in two and three dimensions. In d — 2 
we find a second order transition, in contrast to previous 
results for the bond fluctuation model j2^. We do not 
know whether this is due to an inherent difference in the 
models. In d = 3 a clear first-order phase transition is 
observed, in agreement with all previous studies. But it 
seems that the present study is the first with a detailed 
study of the transition region and of finite size effects. 

We indeed found the latter to be extremely important, 
in particular for d = 3. Estimates of the critical force 
obtained without careful extrapolation to iV — > cx) would 
be grossly wrong. This is similar to the O collapse of un- 
stretched polymers in high dimensions [ p5| where these 
finite size effects even mimicked a first order transition, 
while the true transition is second order. We believe that 
we can exclude the latter for the present case. 

The most precise estimates of critical forces were ob- 
tained via a direct comparison of the free energies in the 
two phases. This is easily done with PERM, in contrast 
to most other Monte Carlo methods. We foimd that the 
partition sum of the collapsed phase is very closely re- 
lated to that of completely unstretched collapsed poly- 
mers. An essential part of our numerical effort went in- 
deed into improved estimates of the latter. In particular, 
we verified that the concept of a surface tension applies 
both to 2-d and 3-d collapsed polymers, and for d = 3 
we verified its scaling with temperature. 

We compared the finite N corrections in detail with 
phenomenological models which yield different order 
transitions in d = 2 and d = 3, and found perfect agree- 
ment. The latter is a bit surprising since fluctuations 
and surface layer thicknesses are neglected in the latter. 
These approximations should become exact in the limit 
of chain length iV — )■ cxd, but our present chains, with 
N w 10"^ — 10'', could be expected a priori to be much 
too short for this. 



APPENDIX 

PERM is a particular depth-first implementation of 
sequential importance sampling with reweighting [^ . 
Polymer chains are built like random walks by adding 
one monomer at each step. 

As in any such algorithm we have the freedom to sam- 
ple these steps from a wide range of possible distribu- 
tions, provided this additional bias is taken into account 
by suitable weight factors. First of all, we used a Rosen- 
bluth like bias in avoiding steps which would lead to 
self intersections. Among the other possible steps we 



selected those parallel, antiparallel, and transverse to F 
with probabilities p+x ■ P-x '■ P± = Vb : a/1/6 : 1. 
The Boltzmann factors for the pair interactions were 
taken into account entirely by the weight factors and 
did not enter into the step probabilities. Let us define 
— 1, pf^ — h'^^l'^ . Then we have = if step i is 
forbidden, and 



ii— /allowed J 

else. The corresponding weight factors are then 

w, = (Aa;, = 0, 1, or - 1) , (25) 

Pi 

where r7i„ is the number of neighbours of the new site 
already occupied by non-bonded monomers. 

The total weight of a chain of length n is then Wn — 
Iln'<n^i„'- Every time an 7i-th monomer is added to 
the chain, we update the current estimate of the parti- 
tion sum to 

Z„ = M-i^iy„(a), (26) 

a = l 

where M„ is the number of chains reaching length n and 
Wn [a) is the weight of the ath chain. Chains are cloned 
and pruned if their weight is above iZn and below Zn/i, 
respectively. 

For different sets of simulations we measured different 
observables. In the dense limit simulations, e.g., we mea- 
sured only the partition sum, while the largest number of 
observables was measured for stretched chains on infinite 
lattices. There we measured the partition sum Z„, the 
average end-to-end displacement {x) parallel to the force, 
the average squared end-to-end distance < >, and 
the number of contacts between non-bonded monomers 
< TO > which is a measure for the internal energy. We 
also measured histograms of the parallel displacement 
and of the contact number. In addition during all runs 
we made also some technical control measurements. Most 
importantly, we generated histograms of tour weight dis- 
tributions in order to test whether the results are 
statistically reliable or not. 

In all cases we used a single integer to label lattice sites, 
and used helical boundary conditions. If the lattice size 
was a power of 2^, say 2*^'' with integer k, then the neigh- 
bours of site i are «± 1, i ± 2*^, . . . i ± 2(''~i)'^, and integers 
outside the range [0, . . . 2'^'^ — 1] are brought back to this 
interval by means of bitwise AND with 2'='' - 1. If the 
lattice size was a different power of 2, say 2^, the num- 
bers 1, 2*"', . . . 2('^~^)'^ are replaced by the integers nearest 
to l,2P/'',...2(''-i^P/''. 

In d = 2 and for finite lattices, self avoidance and con- 
tacts were simply checked by means of bit maps: in an 
array of characters, each occupied site was marked '1', 
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while each empty site was marked '0'. For opens systems 
in d = 3 this would have needed too much storage and we 
used hashing. In this way we could implement effectively 
infinite lattices with relatively small computer memory. 
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